Nonlinear Band Structure in Bose Einstein Condensates: The Nonlinear Schrodinger 
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All Bloch states of the mean field of a Bose-Einstein condensate in the presence of a one dimen- 
sional lattice of impurities are presented in closed analytic form. The band structure is investigated 
by analyzing the stationary states of the nonlinear Schrodinger, or Gross-Pitaevskii, equation for 
both repulsive and attractive condensates. The appearance of swallowtails in the bands is examined 
and interpreted in terms of the condensates superfluid properties. The nonlinear stability properties 
of the Bloch states are described and the stable regions of the bands and swallowtails are mapped 
out. We find that the Kronig-Penney potential has the same properties as a sinusoidal potential; 
Bose-Einstein condensates are trapped in sinusoidal optical lattices. The Kronig-Penney potential 
has the advantage of being analytically tractable, unlike the sinusoidal potential, and, therefore, 
serves as a good model for experimental phenomena. 

PACS numbers: 03.75.Hh,05.03. Jp,03.65.Ge 
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I. INTRODUCTION 



Periodic potentials are ubiquitous in physics, appear- 
ing in electron transport in metals 0, Josephson junc- 
tion arrays 0, nonlinear photonic cry" als and waveguide 



rays 



arrays |jg, and Bose-Einstein condensates (BEC's) |^. 
With the realization of BEC's of alkali atoms in a si- 
nusoidal optical lattice, there has been an explosion in 
studies of BEC's in periodic potentials, both experimen- 
tally and theoretically d S S H E3 ■ BEC's in peri- 
odic potentials have been used to study phase coherence 
of atom lasers 0, and matter- wave diffraction 01 • 
Therefore, in the context of the BEC, the study of peri- 
odic potentials provides an excellent connection between 
condensed matter physics and atomic physics. In con- 
trast to other physical contexts, the lattice geom etry and 
strength, as well as the interatomic interactions |14t lla] , 
are all tunable parameters for the BEC. We show that the 
mean field Bloch states of a BEC in a Kronig-Penney po- 
tential, i.e., a lattice of delta functions, exhibits the same 
band structure and stability properties as the experimen- 
tal case of a sinusoidal potential. Unlike in the case of the 
sinusoidal potential, Bloch state solutions to the Kronig- 
Penney potential can be described by straightforward an- 
alytic expressions. The Kronig-Penney potential has dis- 
tinct advantages as a model of BEC's trapped in periodic 
potentials. 

Specifically, we consider the steady state response of 
the mean field of a BEC to a Kronig-Penney poten- 
tial using the Bloch ansatz. The mean field is modeled 
by the nonlinear Schrodinger equation (NLS), which ap- 
pears in numerous areas of physics; in the context of 
the BEC, it is often called the Gross-Pitaevskii equa- 
tion |TfiL . We have previously obtained the full set 
of stationary solutions for a single delta function ana- 
lytically [H H IH EH- Using these results, we are 
able to rigorously describe the band structure. We find 
that above a critical nonlinearity, swallowtails, or loops, 



form in the bands [22, [23 ■ These swallowtails are re- 
lated to superfluid properties of the BEC [24|, as we 
shall explain. Stability properties are studied numeri- 
cally by time evolution of perturbed stationary state so- 
lutions to the NLS. It is found that stable, as well as 
unstable, regimes exist for both repulsive and attractive 
BEC's. We also elucidate the linear Schrodinger and dis- 
crete nonlinear Schrodinger limits of our model. The 
latter limit is related to the superfluid phase of the Bose- 
Hubbard model pa) . 

Experimentally, in order to create a BEC in a lat- 
tice, alkali atoms are first cooled to a quantum degener- 
ate regime by laser co oling an d evaporation in harmonic 
electromagnetic traps |2(| [27| . A sinusoidal optical lat- 
tice is then created by the interference of two counter- 
propagating laser beams, which creates an effective si- 
nusoidal potential proportional to the intensity of the 
beams 10] . The potential is due to the AC Stark shift 
induced by the dipole interaction with the electromag- 
netic field on the atoms' center of mass [28|. For large 
detuning of the optical field from the atomic transition, 
dissipative processes, such as spontaneous emission, can 
be minimized and the potential becomes conservative. 
Nonzero quasimomcntum can be examined by slightly de- 
tuning the two lasers by a frequency 8v [|| . The resulting 
interference pattern is then a traveling wave moving at 
the velocity, v = {\/2)8v, where A is the wavelength of 
the first beam. This produces a system with quasimo- 
mentum q = mv/h, where m is the atomic mass. After 
a given evolution time, the traps are switched off. The 
BEC is allowed to expand and the density is then imaged. 

The sinusoidal optical lattice potential is composed of 
a single Fourier component. If more counter-propagating 
laser beams of different frequencies are added, more 
Fourier components are introduced, and the potential be- 
comes a lattice of well separated peaks. In the limit that 
the width of the these peaks becomes much smaller than 
the healing length of the BEC, the potential effectively 
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becomes a Kronig-Pcnney lattice. 

We proceed to highlight a few of the many experiments 
on BEC's trapped in lattice potentials. In a work by 
Anderson and Kasevich ^lj , the tunneling of a BEC be- 
tween sites of an optical lattice aligned with the gravi- 
tational field was examined. This produced the matter 
wave analog of the AC Josephson effect, as well as a new 
kind of atom laser. Greiner et al. |ld| observed a quan- 
tum phase transition in a 3D lattice from a superfluid to a 
Mott-insulator. BEC's trapped in lattices have been pro- 
posed as one possible realization of a quantum computer 
and many studies have been performed in this direction, 
see for example pflj and references therein. An optical 
lattice was used to renormalize the effective mass to cre- 
ate a Tonks-Girardeau gas [30l l3lL l32| , which describes 
a truly one dimensional Bose gas. With regards to band 
structure, Fallani et al. investigated the dynamic in- 
stability of a BEC in a lattice for various quasimomcnta. 
By determining the rate of loss due to heating at a given 
quasimomentum, they were able to compare their data 
with theoretically predicted instability growth rates. We 
emphasize that our study treats the quasi one dimen- 
sional regime in which the BEC is always in a superfluid 
state and can be modeled by the NLS with a periodic 
potential, as was the case in the experiment of Fallani et 
al. 



The case of a BEC trapped in a sinusoidal potential 
has been studied theoretically in grea t detail by a num- 
ber of researchersJjjJplJp, |M M HI S HE 0, IS 
S ES IS E ElElEir Anah/tic results were re- 
stricted to certain special cases [4(J, |43, 13, ; the full 
band structure was obtained numerically, as for example, 
by expansion of the wavefunction in Fourier modes [22j. 
Many of these authors considered stability properties of 
stationary states. We will compare our results to those 
of Machholm et al. j^j. A generalized periodic poten- 
tial in the form of a Jacobian elliptic function has been 
studied with mathematical rigor |43l 143, l45j . A subset 
of the solutions to the Kronig-Penney potential has been 
found [46l l47ll48Ll49ll50l |. Recently, for example, Li and 
Smerzi |4g , investigated generalized Bloch states for con- 
stant phase and zero current. 

The article is organized as follows. In Sec.[n] the Bloch 
wave solutions to the stationary NLS with a Kronig- 
Penney potential are presented. The quasimomentum- 
cnergy bands are detailed for repulsive and attractive 
condensates in both the weakly and strongly interacting 
regimes in Sec. IIIII and Sec. lIVI In Sec. El the changes in 
the density profile of the condensate is examined as the 
quasimomentum changes. In Sec. I VII various limiting 
cases of the NLSE are examined. The stable and unsta- 
ble regimes of the bands are studied in detail for both 
repulsive and attractive condensates in Scc lVIII Finally, 
concluding remarks are made in Sec. IVIIII 



II. THE NONLINEAR SCHRODINGER 
EQUATION AND BLOCH WAVES 

In this paper, we consider the mean-field model of a 
quasi-lD BEC in the presence of a Kronig-Penney po- 
tential, 
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V(x) = V J2 S(x-jd) 



(1) 



] = -co 



where d is the lattice spacing and Vq is the strength of the 
potential. When the transverse dimensions of the BEC 
are on the order of its healing length, and its longitudinal 
dimension is much longer than its transverse ones, the 
ID NLS [5l| which describes the stationary states of the 
mean field of a BEC is given by, 
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^ xx + g\^\ 2 ^ + V(x)^ = ^ . 



(2) 



Note that harmonic oscillator confinement in the trans- 
verse directions with frequency u> has been assumed for 
atoms of mass m and u is the eigenvalue, g character- 
izes the short range pairwise interaction, and V(x) is an 
external potential |52j . In this paper, a weakly inter- 
acting system is defined by, gn/Vo < 1, and a strongly 
interacting system is defined by, gn/Vo 3> L In the case 
where the harmonic oscillator length approaches the s- 
wave scattering length, a s , the ID NLS no longer models 
the system and a one-dimensional field theory with the 
appropriate effective coupling constant must be consid- 
ered instead jS] ■ Since a s is on the scale of hundreds of 
angstroms for typical BEC's, this regime is not relevant 
to the present study. 

In Eq. the length is scaled according to the lattice 
spacing, d, and energy has been rescaled by tt 2 /(2Eq), 
where 



En = 



2md? 



(3) 



is the kinetic energy of a particle with wave vector equal 
to that at the boundary of the first Brillouin zone. The 
variables in Eq. @ are defined by, 
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2E 



V(x) 



2E 9 ' 

£rV'(x'/d), 



(4) 
(5) 
(6) 
(7) 



where the primed variables contain the physical units 
of the system. The renormalized ID coupling is, gn = 
2na s w j_hmd 2 / h 2 , where n is the average density per lat- 
tice site. Both attractive and repulsive atomic interac- 
tions, i.e., g > and g < 0, shall be considered. The 
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wave function or order parameter ^>(x 7 t) has the physi- 
cal meaning of, 



^(x, t) = y p(x, t) exp[— ifxt + i(j>(x, t)] 



(8) 



where p(x, t) is the line density and the supcrfluid veloc- 
ity is given by v(x,t) = (Ti/m)d4>(x,t)/dx. 

In addition to the NLS, Eq. J3J), the normalization of 
the wavefunction is given by, 



p(x')dx' , 



(9) 



where n is the number of atoms per lattice site. The 
boundary conditions induced by the Kronig-Penney po- 
tential causes a discontinuity in the derivative of the 
wavefunction across each delta function, 



d x p(j + e) - d x p(j - e) = W Q p(0) , 



(10) 



where j is an integer and e — ► 0. 

A brief review is now given of the general solution to 
Eq. with no external potential. We have previously 
presented a proof that this represents the full set of so- 
lutions for a constant potential [l^. Therefore, by us- 
ing this complete set of stationary state solutions to the 
constant potential case, we can calculate the full set of 
Bloch solutions for a lattice. With the solutions in the 
form of Bloch waves, the relationship between the energy 
per particle and the quasimomentum of the Bloch waves 
is determined. 

The constant potential solutions are the form of the 
density in each lattice site. The density, p, and the phase, 
(p, are, 

k 2 b 2 

p{x) = B H sn 2 (bx + x a ,k) , (11) 
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<f>{x) 



o p{x') 



-da;'. 



(12) 



where sn is a Jacobi elliptic function [53l l54| and the 
density offset, B, the horizontal scaling, b, the transla- 
tional offset, xq, and the elliptic parameter, k, are free 
variables. The Jacobi elliptic functions are generalized 
periodic functions characterized by an additional param- 
eter, k e [0, 1]. In the limit that k — ► and k — ► 1, the 
Jacobi elliptic functions become circular and hyperbolic 
trigonometric functions, respectively. The period of the 
square of the Jacobi elliptic functions is given by 2K(k), 
where K(k) is the complete elliptic integral of the first 
kind Hill. 

Inserting Eqs. ljTT]> and (fT^f) into Eqs. I© and @, with 
V(x) = 0, one finds that the eigenvalue, p, and phase 
prefactor, a, are given by, 



1 



a 



k 2 ) + Wg) , 
B(k 2 b 2 /g + B)(b 2 + Bg). 



(13) 
(14) 



Note that the fact that a only enters into the equations as 



for which a / 0, are doubly degenerate, as ±ct lead to 
the same value of the eigenvalue, p, without otherwise 
changing the form of the density or phase. The boundary 
conditions are used to determine the appropriate values 
of the variables. 

After using the methods from the previous paper |l8j | , 
the solutions are now represented in the usual Bloch 
form, 



(15) 



where q is the quasimomentum and f q (x) has the same 
period as the lattice, f q (x) = f q (x + 1). By substitut- 
ing Eq. ijT5)) into Eq. (JSJ), one finds that the density, p, 
must also have the same period as the lattice and that 
the quasimomentum and energy per particle can be de- 
termined from the density profile by, 



o P(%') 



dx' . 



(16) 
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m* + V 6(x)W (17) 



The quasimomentum is simply the phase jump across 
each lattice site and corresponds to the momentum due 
to the supcrfluid velocity of the system. Since nontrivial 
phase solutions are degenerate for the two phase prcfac- 
tors, ±a, only half of a Brillouin zone needs to be cal- 
culated, i.e. < q < ir. The second half will then be a 
reflection around q = 0. 

This paper examines the quasimomentum-energy 
bands and, therefore, reduces the problem to a situation 
where the density is symmetric around the ends, x = j 
or x = j + 1, and the middle, x = j + 0.5, of the lattice 
sites, where j is an integer. Due to this symmetry, there 
are only two possible values for the translational offset, 
xq, 



x &{- h -,K{k)- h -}, 



(18) 



a implies that all nontrivial phase solutions, 



those 



where K(k) is the half period of the density. The offset 
forces the density in the center of each site to be either 
a minimum or a maximum of the site, depending on the 
sign of the interaction. 

Since it is computationally intensive to include the in- 
tegral in the quasimomentum equation, Eq. 1)16(1 . with 
a root finding algorithm, one of the parameters, b, k, 
or B, is varied while the other two are determined from 
the number equation, Eq. 0, and boundary condition, 
Eq. i[TU)l. The quasimomentum and energy are evaluated 
from these parameters and can then be plotted paramet- 
rically. In the following two sections, we discuss the en- 
ergy bands from repulsive and attractive interatomic in- 
teractions. 



III. REPULSIVE ATOMIC INTERACTIONS 

The structure of the energy bands is strongly depen- 
dent on the strength and sign of the atom-atom inter- 
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FIG. 1: Energy per particle as a function of the quasi- 
momentum for the first three bands of a weakly repulsive 
condensate (gn = Eo) in a repulsive lattice (Vb = Eo). The 
noninteracting linear band structure is given by the dashed 
curves. 



FIG. 3: Energy per particle as a function of the quasi- 
momentum wave number for a noninteracting system with 
no potential (dashed curve), a periodic potential and a small 
or no interaction (dot-dashed curve) and a periodic potential 
with a large interaction (solid curve). 
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FIG. 2: Energy per particle as a function of the quasi- 
momentum for the first three bands of a strongly repulsive 
condensate (gn = 10i?o) in a repulsive lattice (Vb = Eo). 



actions, g. In Fig. 2] and Fig. [21 the energy bands for 
specific cases of weak and strong repulsive interactions 
are presented, gn = Eq and gn = 10.Eo, respectively. 
The condensates are assumed to be in a repulsive lattice, 
Vq = Eo. Note that the energies are given as unsealed 
values so that easy comparison with previous literature 
is possible. 

Notice that in Fig. ^ the interaction strength is small 
and deviations from the linear band structure are small 
as well. The bands are vertically shifted higher as com- 
pared to a linear system due to the repulsive interactions 
that increase the energy of the system. When the inter- 
action strength is further increased the band structure 
becomes quite different. Swallowtails [3^ appear at the 
ends of the bands, as in Fig. [21 The width of these swal- 
lowtails grows as the interaction strength is increased. 
Swallowtails are a general feature of a nonlinear system 



in a periodic potential |42j and appear for both repulsive 
and attractive interactions. 

The presence of these swallowtails is due to the hys- 
teric behavior of the superfluid condensate. A thorough 
discussion of this topic is given by Mueller [24| (see also 
references therein). For a completely free, noninteract- 
ing system, the energy has a quadratic dependence on 
the momentum, shown as the dashed parabolic curves 
in Fig. Since a quasimomentum of 2-7T leaves the sys- 
tem unchanged, the quadratic dependence is repeated 
centered around integer multiples of 2tt. When a peri- 
odic potential is added to the system, bands in the en- 
ergy are formed, shown as the dot-dash sinusoidal curve 
in Fig. |3J These bands may be found by solving the 
linear Schrodinger equation. An interacting condensate, 
however, is a superfluid and can therefore screen out the 
periodic potential |4lj . The energy band will then ap- 
pear similar to that of a free particle, shown as a solid 
swallowtail curve in Fig. [31 until a critical point. Since 
there is a critical velocity, determined by the condensate 
sound speed, the energy band must terminate. If this 
velocity allows the quasimomentum to pass the edge of 
the Brillouin zone, there are then two separate energy 
minima. This demands that there be a saddle point sep- 
arating them and hence the three stationary states. For 
a Kronig-Penney potential, just as for a sinusoidal poten- 
tial, when a swallowtail appears there are two spinodal 
points, where the number of energy extrema changes, at 
the edges of the swallowtail. These are points where the 
local minimum and local maximum energy converge to 
the same energy. As the interaction strength is increased, 
the width of the swallowtail can increase without limit. 
Eventually, the swallowtail will become large enough to 
cross the higher bands. In this case, there is a degeneracy 
between bands. 

For interacting systems in periodic potentials, there is 



a minimum interaction strength for which the swallow- 
tails in the energy bands can exist [24[. This can, in 
general, be dependent on both the strength of the poten- 
tial as well as the band that is being discussed. For an 
optical lattice it was shown that the onset of the swallow- 
tail for the lowest band of a repulsive condensate occurs 
when the interaction strength and the potential are equal. 
For hi ghe r bands the relationship no longer becomes an- 
alytic I23. For the Kronig-Penney lattice potential, the 
critical value for the onset of the swallowtails is not de- 
pendent on the band under consideration or the sign of 
the interaction. Numerically, we are able to determine 
that the onset occurs when gn = 2Vq. The factor of two 
is not present for the optical lattice potential. This holds 
true in all situations except for the lowest band of an at- 
tractive condensate, as will be explained in the following 
section. 

The energy bands arc slightly different when the con- 
densate is in an attractive potential, i.e. Vq < 0. At the 
Brillouin zone boundary, the energy gap between bands is 
proportional to Vq for a weakly interacting system. The 
density is an equally spaced array of solitons. With a 
quasimomentum of q = tt in the lowest band, the conden- 
sate has one dark solitonj^jj per lattice site with nodes 
at the delta functions. The potential, therefore, does not 
contribute to the total energy of the system since the 
delta function potential only depends on the density at 
that position. The second band with quasimomentum 
of q = tt, however, has a density that is offset by half 
a period. There is then a dark soliton in the center of 
each lattice site. The effect of the potential lowers the 
energy by p(Q)Vq. The second band then can cross the 
first band for a condensate in an attractive potential. In 
contrast, the bands are separated by an energy p(0)Va 
for a repulsive potential. As the interaction strength in- 
creases, the effects of the potential become less notice- 
able and the bands are no longer degenerate. Note that 
repulsive and attractive sinusoidal potentials create the 
same band structure. The difference in the current sys- 
tem arises since there are two length scales associated 
with a Kronig-Penney potential, the lattice spacing and 
delta function width. 



IV. ATTRACTIVE ATOMIC INTERACTIONS 

The energy bands for an attractive condensate in a re- 
pulsive potential with a small interaction strength have a 
qualitatively similar form as for a weakly repulsive con- 
densate. Note that the attractive bands are, however, 
lower in energy than the repulsive bands. The attrac- 
tive bands are pushed down from the linear case due to 
the attractive interaction strength. A strongly attractive 
condensate, however, has several qualitative differences 
compared to a strongly repulsive condensate. 

In Fig. 01 the band structure for a small attractive in- 
teraction, gn = —i?o, in a repulsive potential, Vq = Eq, 
is illustrated. The band structure for this system is al- 
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FIG. 4: Energy per particle as a function of the quasi- 
momentum wave number for the first three bands of a weakly 
attractive condensate [gn — —Eo) in a repulsive lattice 
(Vb = Eo). The noninteracting linear band structure is given 
by the dashed curves. 



most identical to that of the weakly repulsive condensate, 
Fig. CI This is to be expected since a weak interaction 
should only cause small perturbations to a noninteracting 
system. 

In Fig. [SJ the band structure for a large attractive in- 
teraction, gn = — IOEq, in a repulsive potential, Vq = Eq, 
is illustrated. The energy band for this system is vastly 
different from the strongly repulsive case. The swallow- 
tails in the bands are now on the upper band at the band 
gaps as opposed to the lower bands at the band gaps as 
they were for the repulsive case in Fig. El The first energy 
band never has a swallowtail. This is because the swal- 
lowtail must be on the lower portion of the band and 
below the center of the band there is no quadratic en- 
ergy dependence for the swallowtail to follow, see Fig. 
Therefore, the lowest energy band only has a swallowtail 
for a condensate with repulsive interactions. 

The second band of a strongly interacting attractive 
condensate looks quite different than that of a strongly 
interacting repulsive condensate. For an attractive con- 
densate, after the initial critical value of the interaction 
strength is reached, a swallowtail in the band starts to 
form, as in the third band of Fig. [5] As the interac- 
tion strength increases, the width of the swallowtail also 
increases. Eventually, another critical value is reached 
where the width of the swallowtail in q reaches 7r and 
runs into the band edge. If the interaction strength is 
further increased, the width of the swallowtail increases 
but the quasimomentum that should be less than zero 
becomes imaginary due to the form of a. This represents 
a nonphysical solution and is contrary to the assump- 
tion that the phase is real in Eq. |JS}. Therefore, part of 
the band is absent and the band appears as two sepa- 
rate pieces, a loop and a separate line. These are both 
marked as Band 2 in Fig. [SJ The band takes on this ap- 
pearance since the swallowtail is attempting to extend 
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FIG. 5: Energy per particle as a function of the quasi- 
momentum wave number for the first three bands of a strongly 
attractive condensate (gn = — 10.Eo) in a repulsive lattice 
(Vo = Eo). The first band does not have a swallowtail. How- 
ever, the second and third bands do have swallowtails. The 
swallowtail in the second band appears as a loop with an 
unattached curve since an attractive condensate has a max- 
imum width for the swallowtails, which in the case of the 
second band is a width of tt. 



lower than the minimum in the adjacent quadratic en- 
ergy dependence. In a similar fashion to the first band, 
the solutions are not physical. The portion of the swal- 
lowtail that is before this edge is physical, however, and 
so a third energy extremum is still present. The third 
band will exhibit the same behavior when the width of 
the swallowtail reaches 27r. In general, the n th band will 
appear similar when the swallowtail reaches a width of 
(n — l)7r since the swallowtail will then have reached the 
minimum of the corresponding free particle energy. It 
should be noted that this phenomenon does not occur 
for a repulsive interaction since there is no extremum to 
limit the growth of the swallowtail. 

The energy bands arc slightly different when the con- 
densate is in an attractive potential. In a manner similar 
to a repulsive condensate in an attractive potential, for a 
sufficiently weak interaction strength and large potential 
strength, the first and second bands can become degen- 
erate. 

It should be noted that for a single attractive delta 
function on the infinite line, there is a critical repulsive 
interaction above which the impurity can no longer bind 
the condensate 0. However, for a periodic system of 
delta functions, there is no longer a critical value. This 
can best be understood by considering a condensate of a 
finite size bound by a single delta function. As the re- 
pulsive interaction increase's, the condensate spreads out 
until it reaches the next delta function. This continues in- 
definitely as the interaction strength increases. There is, 
therefore, no transition from a bound to unbound state as 
there was with a single impurity. The symmetry break- 
ing nature of a single impurity allows for the transition 
that is not present with the symmetric lattice. 



V. DENSITY PROFILES 

As density is the primary experimental observable for 
BEC's, the change in the density profile as the quasimo- 
mentum is varied is important. The case of weak inter- 
actions creates similar changes in the density profile and 
energy band structure independent of the sign of the in- 
teraction. This is because the interaction energies are of 
the same magnitude as the potential. For strong interac- 
tions, the density is in general more sharply peaked for 
an attractive condensate and flatter for a repulsive con- 
densate. However, the general way in which the density 
changes is qualitatively similar. 

The density changes in the first band of the weakly at- 
tractive condensate are shown for three different quasi- 
momenta in Fig. EJ The solid curve, dashed curve, and 
dotted curve represent the density profile for q = 0, 
q = 7r/2, and q = tt, respectively. To understand the 
energy bands in terms of the density profile, the three 
terms of the energy in Eq. ifTTjl should be discussed. The 
kinetic, interaction and potential energy per particle are 
given by, 
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|l*| 4 dZ: 



El = I f V 6(x)\#\ 2 dx . 
n n J 



(19) 
(20) 
(21) 



respectively. Notice that the density at the origin mono- 
tonically decreases as the quasimomentum is increased. 
Therefore the potential energy will also monotonically 
decrease due to the delta function at x — j, where j is 
an integer. Since the condensate is attractive, the inter- 
action energy decreases monotonically as the quasimo- 
mentum is increased since the density is becoming more 
peaked. The kinetic energy monotonically increases as 
the quasimomentum increases since the variations in the 
density become larger. For a weakly repulsive conden- 
sate, the density will also become more strongly peaked 
as the quasimomentum is increased and will thus have 
an interaction energy that increases. The larger varia- 
tions in the density will then increase the kinetic energy 
as well. 

The case of the density variations associated with the 
first band of a strongly interacting repulsive condensate 
is shown in Fig. [7] The change in the density is qualita- 
tively similar to that of the weakly attractive case. The 
solid line represents the density profile for q = 0. The 
dot-dash curve represents the density when q = tt at the 
bottom of the swallowtail. The dashed curve represents 
the density when q = 0.467T at the end of the swallow- 
tail. The dotted curve represents the density when q = tt 
at the top of the swallowtail. The density at the ori- 
gin again monotonically decreases as the quasimomen- 
tum increases and therefore the potential energy also in- 
creases monotonically. The interaction energy decreases 
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FIG. 6: Changes in the density of a weakly attractive, gn — 
—Eo, condensate associated with different positions on the 
first band. The solid curve represents the density when q = 0. 
The dashed curve represents the density when q = n/2. The 
dotted curve represents the density when q — it. The lower 
plot shows the corresponding positions on the first energy 
band. 



monotonically as the band is traversed since the density 
becomes more uniform. The kinetic energy follows the 
same qualitative path as the total energy. It is therefore 
the kinetic energy that has the greatest influence on the 
energy bands. 

In the case of strongly attractive interactions, the low- 
est band consists of well separated bright solitons cen- 
tered in each lattice site. As the quasimomcntum changes 
from zero to 7t, the density changes very little. This is 
easily noted in the energy spectrum of the lowest band 
which changes on the order of a tenth of a percent [3] 



VI. LIMITING CASES 

It is now shown how the solutions to the NLS con- 
nect to the solutions of the linear Schrodinger equation 
as well as to the discrete nonlinear Schrodinger equation. 
The linear Schrodinger equation is important since it de- 
scribes electron motion through crystals. The discrete 
nonlinear Schrodinger equation describes superfluid sys- 
tems where the potential strength is much greater than 
the interaction strength. It is a limit of the Bosc-Hubbard 
model when the Hamiltonian is projected onto a coherent 
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FIG. 7: Changes in the density of a strongly repulsive con- 
densate, gn = WEo, associated with different positions on the 
first band. The solid curve represents the density when q — 0. 
The dot-dashed curve represents the density when q = n at 
the bottom of the swallowtail. The dashed curve represents 
the density when q = 0.46-7T at the end of the swallowtail. The 
dotted curve represents the density when q = n at the top of 
the swallowtail. The lower plot shows the corresponding po- 
sitions on the first energy band. 



state 1 



A. Linear Schrodinger Equation Limit 

The energy dispersion relation of the noninteracting 
regime can be determined from the solutions of the in- 
teracting system. In the linear limit, the interaction 
strength, g, and elliptic parameter, k, must vanish such 
that, 



lim 



k 2 b 2 
9 



A. 



(22) 



where A is the height of the density fluctuations. The 
elliptic parameter must vanish since in this limit the 
Jacobi elliptic functions become circular trigonometric 
functions, which is the appropriate linear form of the 
density. The linear solutions are then of the form, 



p(x) = B + Asm 2 (bx + jir /2 - 6/2) , 



(23) 



where j is a positive integer and b is the horizontal scal- 
ing. The offset jw/2 can be dropped, since the effect in 
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the square of the sine function is to change it into ei- 
ther the square of a sine or a square of a cosine and this 
merely alters the values of A and B. The quasimomen- 
tum equation, Eq. (|16[) . can be integrated exactly in the 
linear limit to find 



cos 2 (<?/2) = 



1 



1 + (1 + r) tan 2 (6/2) ' 



(24) 



where r = A/B is the ratio of the density modulation 
coefficient, A, to the constant density offset, B. The 
equation that governs the discontinuity in the derivative 
of the density at the delta functions, the boundary con- 
ditions given by Eq. I|1(J|) . can be solved for the ratio r, 



V 



6sin(6/2) cosO/2) + V sin 2 (6/2) 



(25) 



The quasimomcntum can then be related to the horizon- 
tal scaling, b, through the relation, 



V 

cos(g) = cos(fr) + — sin(6) . 

b 



(26) 



This dispersion relationship coincides with that of the 
linear Schrodinger equation in introductory quantum me- 
chanics texts |56|. At the edges of the Brillouin zone, 
q = and q = tt, the heights of the bands and the gaps 
can be extracted. 



B. Discrete Nonlinear Schrodinger Equation Limit 

It is now shown how the solutions to the NLS con- 
nect to the solutions of the discrete nonlinear Schrodinger 
equation (DNLS). The DNLS describes a system where 
the density is localized in predominantly non-overlapping 
regions described by the Wannier functions Wj, 



Wj{x-Xj) = —== 



"£(x)dg, 



(27) 



where Xj is the position of the j lattice site, N is the to- 
tal number of lattice sites, and the integral is over the first 
Brillouin zone. Since the regions are non-overlapping, it 
is required that, 



Wj(x)wk(x)dx w Sjk 



(28) 



where Wj(x) is the Wannier function localized at site j 
and Sjk is the Kronecker delta function Q. This is also 
known as the tight-binding approximation |57| . In addi- 
tion, it is assumed that particles can only hop to a nearest 
neighbor site, i.e. particles in the j th site can only move 
to the (j ± l) th site. Therefore, all integrals of the form, 



(d x w*)(d x w k )dx , 



(29) 



vanish unless fc £ {j,j ± I}. One possible way to create 
a system like this is to make the strength of the poten- 
tial sufficiently strong in comparision to the interaction 
strength of the condensate, Vo/gn 1. 

The DNLS may be extracted from the NLS by using a 
single band approximation and writing the wavefunction 
as a linear combination of site-localized Wannier func- 
tions, 



y(x,t)=J2'<Pj(t)Wj(x): 



(30) 



where tpj (t) is the probability amplitude to find a particle 
at site j and Wj (x) is the density profile that is localized 
on site j. When Eq. l|5U|) is substituted into Eq. and 
the spatial degrees of freedom are integrated over, the 
evolution of the system becomes governed by, 

idtijjj = ejijjj + Jjtyj+i + tpj-i) + Ujl^j] 2 ^ , (31) 

where ej is the on-site energy of the j th site, J is the 
hopping strength, and U is the interaction strength, with 
each defined by, 



Vo\w j (0)\ 2 + ± j \d x w,\ 2 dx. 
(d x w*)(d x Wj±i)dx , 



Uj =9 \wj\ dx 



(32) 
(33) 
(34) 



The stationary states of the DNLS are given by plane 
waves, tpj = ipoe 1 ^^^, where p is the momentum and 
the eigenvalue, fi, is given by, 



pL = U\ip \ 2 + e + 2,7 cos(p) . 



(35) 



The subscript j has been dropped since all Wannier func- 
tions have the same profile, only translationally shifted 
to the appropriate site. The energy of the condensate for 
each site is given by, 



|Vo| 2 e + 2|Vo| 2 Jcos(p), 



(36) 



If the Wannier functions are normalized to unity, then 
ipo = \pn- and the equation for the energy per particle 
becomes, 



E s 



1, 



-UN 



2 Jcos(p) 



(37) 



The Wannier functions for the Kronig-Penney potential 
can be extracted from Eq. (|ll|l by taking the limit when 
the potential strength approaches infinity. 

It should be noted that a large attractive interaction in 
a repulsive potential will have a similar effect as a large 
repulsive potential and hence the DNLS can be used to 
describe the system. In this case, the large attractive 
interaction causes the density of each site to approach a 
hyperbolic secant squared profile, i.e., a bright soliton. 
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Hence, k — > 1 and B — > —b 2 /g. The normalization of 
the density fixes 6 = gn/2 and the quasimomentum be- 
comes zero. This density profile is the equivalent of the 
Wannicr functions discussed earlier. The on-site energy 
is calculated to be, 



e = —gnVo exp[gn/2] 
The interaction term becomes, 



2 2 

g z n z 
24 



U 



1 2 

g z n z 
6 



(38) 



(39) 



With the assumption that |g| is large, the hopping term 
is given by, 



J 



9 2 n 2 



(4 + gn) exp[gn/2] 



(40) 



The DNLS energy, Eq. , becomes consistent with the 
NLS energy, Eq. I|17(l . to within three tenths of a per- 
cent for the first band of the strongly attractive case, 
gn = — IOE and Vb = Eq. When the on-site energy, in- 
teraction energy, and hopping strength are all assumed to 
remain constant, the maximum error in energy remains 
within three tenths of a percent across the entire band. 
The energy is slightly underestimated at zero quasimo- 
mentum and is slightly overestimated at 7r quasimomen- 
tum. In this instance, knowing the density profile at one 
value of the quasimomentum allows for extension to the 
full spectrum of quasimomentum. When the interaction 
energy becomes strong enough that only terms to high- 
est order of g need to be considered, the total energy 
becomes —g 2 n 2 /2A and the variation across the band be- 
comes — g 3 n 3 /2 exp[<?n/2], which is essentially flat. 



VII. STABILITY 

We proceed to study the stability of the Bloch states 
and to determine the stable regions of the bands are 
determined. In addition to stable solutions, solutions 
that have instability times much longer than experimen- 
tal time scales can be observed in experiments. Recent 
studies of the stability of condensates in a periodic po- 
tential have focused on linear energetic and dynamic sta- 
bility, also called Landau stability H IH |H |H|5|| . In 
contrast, we consider the full response of the conden- 
sate to stochastic perturbations. In order to numerically 
simulate the NLS with a periodic potential, a ring ge- 
ometry is used with a quantized phase. To ensure that 
the phase quantization does not effect the stability prop- 
erties, enough lattice sites were used to allow for many 
rotations of the phase such that Nq = 2nj, where q is 
the quasimomentum, N is the number of sites nad j is 
an interger. The outcome of the stability analysis is in- 
dependent of the number of sites for sufficiently large 
number of sites. In most cases, j ' = 4 was found to be 
adequate to extract the correct stability properties. 



The delta functions are simulated by single point dis- 
tortions in the potential grid. They are also implemented 
by using boxes of different widths with their area normal- 
ized to create the appropriate potential strength. The 
size of the boxes did not influence the stability properties 
until the width became approximately 10% of the size of 
the healing length, £ = Ti/^2gn. The NLS is evolved us- 
ing a variable step fourth-order Runge-Kutta algorithm 
in time and a filtered pseudo-spectral method in space. 
The noise introduced into the simulations comes from 
the round off error associated with numerical simulations. 
To ensure that the form of the noise from round off er- 
ror did not affect the stability properties of the system, 
initial stochastic white noise of various levels was intro- 
duced into the Fourier spectrum. For levels significantly 
greater than the round off noise, the stability times ap- 
proached those from the round off noise. Introduction of 
white noise at the level in the eighth significant digit pro- 
duced the same instability times as the round off noise, 
which effects the sixteenth significant digit. All simula- 
tions were performed over time scales longer than exper- 
imental lifetimes of the BEC, which are on the order of 
seconds. 

The time at which the onset of instability occurs is 
determined by the effective variance in the Fourier spec- 
trum, 



a(t) 



2E(/(P,0)) 2 



(41) 



where /(p, t) is the Fourier component of the wavefunc- 
tion at momentum p and time t and the sum is over 
the momentum grid. This quantity determines how dif- 
ferent the Fourier spectrum is compared to the original 
stationary state. It vanishes when the two spectrums are 
identical and approaches unity when there are no Fourier 
components in common. When a(t) reaches 0.5, i.e. 50% 
of the Fourier spectrum is different than the original, the 
system is considered to have become unstable. 

Unless otherwise noted, for the stability analysis the 
lattice spacing is given by d = 1 /zm, the length scale with 
which current optical lattices are created. In addition, all 
instability time scales will be given for 87 Rb. 



A. Attractive Atomic Interactions 

With an attractive interaction of gn — ~Eq in a repul- 
sive potential of Vb = Eq (see Fig. QJ, the lowest energy 
solution, zero quasimomentum in the lowest band, has 
a lifetime greater than experimental time scales. How- 
ever, when even a slight harmonic perturbation to the 
potential is added to the initial time step, for instance 
a harmonic frequency of 120 Hz which is approximately 
the experimental trapping frequency [59j | in the longitudi- 
nal dimension, the condensate becomes unstable on time 
scales on the order of 1.5 ms. This is short compared to 
the lifetime of a BEC [53, but is still observable. 
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When the quasimomentum of the first band is in- 
creased, the stability of the system becomes dependent 
on the effective mass. The effective mass, m*, is defined 
by0, 



(42) 



d 2 E/dq 2 



The effective mass is the mass that the particle appears to 
have if the potential was not being considered The 
sign of the effective mass can be transferred to the in- 
teraction strength, changing an attractive interaction to 
a repulsive interaction. Therefore, when the quasimo- 
mentum increases and the energy band becomes concave 
down, m* < 0, the system enters a regime of stability. 
The system remains stable even in the presence of the 
harmonic perturbation. 

For zero quasimomentum in the second band, the 
system immediately possesses periodic variations in the 
phase and density. There is an additional underlying in- 
stability that occurs on the order of 5 ms that destroys 
the periodicity of the system. It would be expected that 
this part of the band be stable since there is a nega- 
tive effective mass but the oscillations due to the two 
bright solitons per lattice site force the system to be un- 
stable. The oscillations, although periodic in time, create 
a larger underlying instability to grow. 

The stability properties of a strongly attractive con- 
densate are similar to those of a weakly attractive con- 
densate. The stability of the first band is determined by 
the effective mass while higher bands always go unsta- 
ble. Therefore, for an attractive condensate, the system 
of Bloch waves is stable only if there is one soliton per 
lattice site and the effective mass is negative. 
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FIG. 8: Logarithm of the Fourier spectrum during the time 
evolution. Time is in units of Ti/Eq and distance in units of 
the lattice spacing. 
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B. Repulsive Atomic Interactions 

Similar to an attractive condensate, a repulsive con- 
densate only has stable regions on the first band. For 
a weakly interacting repulsive condensate, gn = Eq and 
Vq = Eo, the effective mass in the first band is positive 
between q — and q = ir/2. The effective mass then 
becomes negative for larger quasimomentum since the 
energy becomes concave down and hence the system be- 
comes unstable. For a quasimomentum of q = 97r/16 the 
instability time is 10 ms and reduces to 2 ms for q = n. 
In this regime, with negative effective mass, the ground 
state is an envelope soliton that can spread over many 
lattice sites. These types of states are called gap soli- 
tons 0, IU 113, E01 an d only occur in interacting systems. 
Figure [5] presents the unstable evolution of the weakly 
repulsive condensate in the first band with a quasimo- 
mentum of q = 7r in Fourier space. Notice that the in- 
stabilities arise from perturbations around the primary 
Fourier components of the wavefunction. In Fig. the 
effective variance, cr, is plotted as a function of evolution 
time and the system becomes unstable around 2 ms. The 
second band becomes unstable when q = in 0.500 ms 



12 3 

t (units of md 2 /R) 

FIG. 9: The time evolution of the effective variance of the 
momentum density, a. Time is in units of h/Eo- 



and in 6 ms for q = tt. Therefore, the system is stable in 
the first band with positive effective mass, m* > 0, and 
unstable elsewhere. This is consistent with the effects 
that the effective mass has on stability in systems de- 
scribed by the lowest band DNLS. In a work by Fallani 
et al. |(|, the instability time of a condensate in a lat- 
tice was measured by using an RF-shield to remove the 
hottest atoms produced by the heating created in the 
sample due to instability. The loss rate, given by the in- 
verse of the lifetime, should then be qualitatively similar 
to the instability time. Our calculations are consistent 
with these experimentally observed loss rates of a BEC 
in an optical lattice 0. 

Due to the presence of the swallowtails, the strongly in- 
teracting system provides different stability regimes. For 
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2 3 
Lattice Spacing (units of i;) 

FIG. 10: The time to instability as the lattice spacing is var- 
ied. Time is in unit of Ti/gn and distance in units of the 
healing length, £ = h/^/2gn. 



a repulsive condensate, gn = IOEq, in a repulsive lattice, 
Vo = E , the main section of the first band, as well as 
the lower portion of the swallowtail, have positive effec- 
tive mass and remain stable. The upper portion of the 
swallowtail, as discussed in Sec. IIIII is an energy maxi- 
mum and is not expected to remain stable. In fact, the 
instability time of the upper portion of the swallowtail is 
approximately 0.2 ms, independent of the actual quasi- 
momcntum. 

The size of the lattice spacing can influence the time 
scales for which the system becomes unstable. In Fig. 1101 
the time to instability is presented as a function of the 
lattice spacing for a repulsive condensate, gn — E a , in 
a repulsive lattice, Vq = E a , with a quasimomcntum of 
q = 7r, where E a = h 2 ir 2 /2m(lfj,m) 2 . 

There is a minimum instability time as the lattice spac- 
ing varies that occurs when the lattice spacing is approx- 
imately twice the healing length, £. The time to insta- 
bility is given by half of the interaction strength time, 
gn/h. When the lattice spacing is much larger than the 
healing length, the density becomes extremely flat except 
at the delta functions, where dark solitons form. Since 
the lattice spacing is large, the dark solitons are far apart 
and are effectively noninteracting pinned solitons. Dark 
solitons are known to be robustly stable 0, |5f| and, 
therefore, for a lattice spacing much larger than the heal- 
ing length the system should become stable. For lattice 
spacing smaller than twice the healing length, the con- 
densate has difficulty distinguishing between the separate 
delta functions and sees closer to a constant potential. 
In this regime, the kinetic energy becomes much greater 
than the interaction energy and potential energy since 
variations in the density occur on the length scale of d/2, 
which is less than the healing length. The system then 
becomes effectively free and noninteracting and, there- 
fore, approaches stability. 



VIII. DISCUSSION AND CONCLUSIONS 

The full set of stationary states of a Bosc-Einstcin con- 
densate in a Kronig-Penney lattice potential, with period 
commensurate with the lattice, have been presented an- 
alytically in the form of Bloch waves for both repulsive 
and attractive interactions. The quasimomcntum energy 
bands were found to exhibit a cusp at the critical interac- 
tion strength gn = 2Vo, where g is the interatomic inter- 
action strength, n is the linear density, and Vq is the lat- 
tice potential strength. For larger interaction strengths, 
swallowtails form in the bands. These swallowtails have 
the same qualitative form as for a sinusoidal potentail 
and exhibit the same stability properties. 

Both attractive and repulsive condensates were found 
to be dynamically stable only in the first band and when 
the effective interaction, sgn(m*)<7, was positive. Also, 
even in the first band, the upper edges of swallowtails 
were always unstable. It should be noted that the effec- 
tive mass docs not influence the stability of the higher 
bands since they are always unstable. Therefore, for an 
attractive condensate, the only stable Bloch states exist 
in the first band between the quasimomcntum where m* 
becomes negative and q = ir. A repulsive condensate is 
only stable in the first band from q = to the quasimo- 
mentum where m* becomes negative. Higher bands will 
always be unstable, for both attractive and repulsive con- 
densates. When solutions became unstable, our numer- 
ically studies consistently observed that the instabilities 
originated around the primary Fourier components of the 
wavefunction. This is in agreement with the formal proof 
of the instability of constant phase Bloch-type solutions 
by Bronski et al. |44|. The instability time was found to 
be a function of the lattice constant. If the delta func- 
tions are spaced either smaller than or much larger than 
the healing length of the condensate, the solutions had 
instability times longer than lifetime of the BEC. Thus 
experiments could access formally unstable sections of 
the energy bands, and, by controlling the ratio of the 
healing length to the lattice constant, directly observe 
the dynamics of instability. The results of our stability 
analysis arc consistent with the experimental work per- 
formed by Fallani et al. [jg, in which the loss rate, the 
inverse of the lifetime, was determined by removing the 
hottest atoms with an RF-shield. 

An interesting phenomenon to note is that for a re- 
pulsive condensate when a swallowtail is present, the 
entire energy band is concave up and, hence, the effec- 
tive mass is always positive. This is in contrast to a 
weakly repulsive condensate, when the concavity of the 
energy band changes, creating a region of negative effec- 
tive mass. Therefore, there is a maximum interaction 
strength for which gap solitons 0, 0, IH, H3, ED can 
be formed since they require a negative effective mass, 
which does not occur when a swallowtail is present. The 
maximum interaction energy is given by gn = 2Vq, the 
strength at which swallowtails appear. 

Stationary solutions to the NLS with a Kronig-Penney 
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potential need not take the form of Bloch states. Solu- 
tions with a period which is an integer multiple of the 
lattice period have been shown to exist for the sinusoidal 
potential [5^ and are expected also to be present for 
the Kronig-Penney potential. Envelope solutions, such 
as gap solitons, also play an important role in other sys- 
tems modelled by the NLS and have been observed in 
BEC's 0- The analytic methods which we have de- 
scribed here are equally applicable to these solution types 
and will form the subject of future study |62l |. 
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